%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% Fonction de calcul de la matrice ie, source du vecteur d'efforts 
% g�n�ralis�s �l�mentaire pour l'estimation par patch flux-free.
%
% Entr�es :
%
% a,b,c,d : [scalaires] coefficients.
% epsilon : =1 sur le premier �l�ment du patch,
%           =-1 sur le deuxi�me �l�ment du patch,
%
% Sorties :
%
% ie : [matrice 4x4] matrice ie

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function [ ie ] = Ie( a,b,c,d,epsilon,poutre,sousPoutre )

% Calcul de iA
%%%%%%%%%%%%%%

iA = zeros(4,4) ;

iA(1,1) = (3/2)*(-b*c+a*d) ;
iA(1,2) = 3*b*c + (3/2)*(a-d)*(c+d) ;
iA(1,3) = (3/2)*(b*c-a*d) ;
iA(1,4) = (3/2)*(a+b)*(c+d) ;

for i=2:2:4
    
    delta=i-3;
    
    iA(i,1) = (1/2)*(-b*(3*c+delta)+3*a*d) ;
    iA(i,2) = b*(3*c+delta) + (1/2)*(a-b)*(3*c+3*d+delta) ;
    iA(i,3) = (1/2)*(3*b*c-3*a*d+b*delta) ;
    iA(i,4) = (1/2)*(a+b)*(3*c+3*d+delta) ;
    
end

for k=1:4
    iA(3,k) = -iA(1,k) ;
end




% Calcul de iB
%%%%%%%%%%%%%%

iB = zeros(4,4) ;

iB(1,1) = 3*epsilon*c*(c-(2/5)*d) ;
iB(1,2) = c*epsilon*(c-(1/5)*d) ;
iB(1,3) = 3*epsilon*c*(c-(3/5)*d) ;
iB(1,4) = -c*epsilon*(c+(1/5)*d) ;

for i=2:2:4
    
    delta=i-3;
    
    iB(i,1) = (1/2)*(-(3*c+delta)+3*a*d) ;
    iB(i,2) = (1/2)*( (3*c+delta)*(a+b) + 3*d*(a-b) ) ;
    iB(i,3) = (1/2)*(3*b*c-3*a*d+b*delta) ;
    iB(i,4) = (1/2)*( (a+b)*(3*c+delta) + 3*d*(a+b) ) ;
    
end

for k=1:4
    iB(3,k) = -iB(1,k) ;
end






% Assemblage
%%%%%%%%%%%%

ie = iA + iB ;

ie=ie/(poutre.L^2*sousPoutre.L) ;

end